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I. OBJECTIVES 


Investigations of the collisional evolution of particulate disks subject to the gravita- 
tional perturbation of a more massive particle orbiting within the disk are underway. 
Both numerical N-body simulations using a novel collision algorithm and analytical 
kinetic theory are being employed to extend our understanding of perturbed disks in 
planetary rings and during the formation of the solar system. Particular problems pro- 
posed for investigation are: (1) The development and testing of general criteria for a 
small moonlet to clear a gap and produce observable morphological features in planetary 
rings; (2) The development of detailed models of collisional damping of the wavy edges 
observed on the Encke division of Saturn’s A ring; and (3) The determination of the 
extent of runaway growth of the few largest planetesimals during the early stages of 
planetary accretion. 


n. PROGRESS TO DATE 
A. NUMERICAL SIMULATION STUDIES 


A computational code is under development to simulate the collisional evolution of 
a particulate disk subject to the gravitational perturbation of a massive body orbiting 
within the disk. To obtain the highest possible resolution, the simulation is confined to a 
narrow ring of particles which is still significantly wider than the radial excursions of the 
particle orbits. The orbits of the perturbing mass and the disk particles are closely spaced 
so that the synodic period between successive encounters with the perturber is many 
times greater than the orbit period. As a result, encounters between disk particles and the 
perturber are relatively short, lasting about one orbital period. For some of our investiga- 
tions, we will adopt Showalter and Bums’ (1982) approach for estimating "instantane- 
ous" changes in the particle orbits due to encounters with the perturber. An important 
new feature in our work, relative to that of previous investigators on this problem, is the 



careful inclusion of particle-particle collisions. We plan to study their important effects. 
The computational code is being developed so that sizable portions of particle orbits are 
advanced in sufficiently short time steps so as to allow the simulation of particle colli- 
sions. 

The N-body computer program under development consists in part of a modified 
version of a computer code developed by Hausman and Roberts (1984a) and Roberts and 
Hausman (1985, 1986) for the study of supersonic rarefied gas flows in gas centrifuges. 
In order to compute particle collisions efficiently, the code adopts Monte Carlo tech- 
niques (Bird 1969, 1976). Trajectories are followed exactly, but not the separation of 
pairs. Instead, at the end of each time step, pairs of particles within each (small) volume 
element (or "cell") are randomly chosen to collide, without regard to their separation. 
The rational for this is the near uniformity of the density and velocity field within a 
sufficiently small volume element, so that each particle represents N real particles spread 
throughout the cell. Variable computational cell volumes are being incorporated into the 
computational code in order to achieve higher resolution in regions where density or 
velocity gradients are large. A numerical scheme for choosing the optimal cellular mesh 
is being developed. The increased efficiency resulting from this effort should permit us 
to simulate 10,000 to 30,000 particles, which is an order of magnitude more particles 
than simulated by previous investigators working on this problem, e.g., Schwarz (1981) 
and Showalter and Bums (1982). 

In addition to the Monte Carlo code described above, we are also adapting portions 
of an N-body program developed by Roberts and Hausman (1984) and Hausman and 
Roberts (1984b) and Roberts and Stewart (1987) for the simulation of gas clouds in the 
interstellar medium of spiral galaxies. This second code has a collision algorithm which 
is slightly more accurate (and time consuming) than the Monte Carlo code since particle 
pairs are required to actually touch in order for a collision to occur. The graphics plot- 
ting capabilities of this code are being refined in order to be particularly useful for 
displaying the results of our simulations of perturbed particle disks. These capabilities 
are expected to provide high resolution and fine detail that can be displayed in our 
results. In particular, well defined shock structures should be able to be resolved in these 
particle simulations. Similarly sharp density contrasts are ubiquitous of perturbed 
regions of Saturn’s rings (e.g. Cuzzi et al. 1984). To fully exploit the advantages of these 
two codes we are currently carrying out work to merge the capabilities of the galaxy 
simulation code with the efficient collision algorithm of the Monte Carlo code. 

For planetesimal disks consisting of kilometer-size particles, gravitational 
encounters among the particles are usually as important as physical collisions (Wetherill 
1980b). However, previous numerical simulations of planetesimal disks by Wetherill 
(1980a) indicate that the small disk particles will generally maintain relative velocities 
comparable to their surface escape velocities so that the two-body approximation for 
mutual gravitational scattering and collision cross sections should be sufficiently accurate 
(Wetherill and Cox 1984). Given our Monte Carlo treatment of particle collisions, we 
believe that sufficient accuracy will be achieved by modeling gravitational scattering as 
perfectly elastic collisions (This is consistent with recent analytical treatments, e.g. Shu 
and Stewart 1985). This, of course, does not apply to the gravitational perturbations due 
to the massive perturber which must be treated by numerical integration of the particle 
trajectories. 



In order to address the problem of runaway growth of planetary embryos, we shall 
also simulate cases where the disk particles are allowed to be accreted by the massive 
perturber. Due to the nature of low-eccentricity three-body orbits, an accurate integra- 
tion of trajectories which collide with the massive body would require a variable integra- 
tion step size (e.g. Weidenschilling and Davis 1983) and possibly excessive computer 
time. In order to avoid this difficulty, we shall explore the possibility of using a statisti- 
cally accurate collision cross section for slow three-body orbits, exploiting as much as 
possible the recent results of Wetherill and Cox (1984, 1985). 


B. ANALYTICAL STUDIES 


Re-examination of Radial Diffusion in Planetary Rings 

According to the kinetic theory of transport processes in planetary rings reviewed 
by Goldreich and Tremaine (1982) and by Stewart et al. (1984) the rate of radial mass 
transport is governed by an effective shear viscosity that is deduced from the off- 
diagonal component of the velocity dispersion tensor. For rings of low optical depth the 
dispersion tensor may be calculated from either the Boltzmann or Krook kinetic equa- 
tions, suitably modified to model the inelasticity of collisions among ring particles (Gol- 
dreich and Tremaine 1978, Shu and Stewart 1985). The theory predicts that planetary 
rings of low optical depth spread radially so as to decrease the average surface density of 
ring particles. The accuracy of this prediction has been called into question by recent 
numerical solutions of the Boltzmann equation that seem to indicate a tendency toward 
radial clumping rather than spreading at low optical depths (Brophy and Esposito 1986). 
The possibility of a diffusional instability at low optical depth is exciting because it 
would provide a much needed mechanism for forming the unexplained structure of 
Saturn’s C ring and the surprisingly sharp density gradients at the inner edges of the A 
and B rings. 

We have begun a re-examination of the kinetic theory of planetary rings in order to 
identify a possible flaw in the published theory and to formulate a new theoretical deriva- 
tion of the radial mass transport rate which avoids such flaws. It has been pointed out by 
Greenberg (1986) that when ring particles with different semi-major axes collide inelasti- 
cally, their center-of-mass velocity is often greater than the circular orbit velocity at the 
smaller semi-major axis or less than the circular orbit velocity at the larger semi-major 
axis. The net result of such collisions is often more closely spaced semi-major axes 
rather than radial spreading. This peculiar reversal of the velocity gradient is a conse- 
quence of the non-local character of eccentric orbits that are associated with the radii of 
their semi-major axes rather than the actual radial coordinate of the particle. Such non- 
local effects may not be adequately modeled by published kinetic theory because it is 
formulated in terms of the local density of particles in phase space. In particular, the 
assumption of a tri-axial Gaussian velocity distribution may preclude the local shear 
reversal described by Greenberg because even functions cannot model a velocity distri- 
bution that is systematically positive when the particle is at radii smaller than the semi- 
major axis and systematically negative when the particle is at radii larger than the semi- 



major axis. 

To avoid this difficulty, we have begun a direct calculation of the radial flux of 
semi-major axes due to particle collisions that does not require the determination of an 
effective shear viscosity from the off-diagonal component of the dispersion tensor. The 
calculation is motivated by analogous calculations of the diffusion of plasma guiding 
centers across magnetic fields (see for example Ichimaru 1973). The basic set-up of the 
calculation is similar to the one described by Hameen-Anttila (1978) where one must find 
the first and second moments of the collisional rate of change of the semi-major axes. By 
relating these quantities to the collisional rate of change of the velocity dispersion tensor, 
Hameen-Anttila obtained an effective shear viscosity identical to that obtained by Shu 
and Stewart (1985). At the present time we are seeking to remove a dubious approxima- 
tion in Hameen-Anttila’ s calculation that is functionally equivalent to the assumption of 
a Gaussian velocity distribution that is made in the usual kinetic theory calculation. Prel- 
iminary results of our investigation suggest that the actual radial mass transport rate may 
differ substantially from the rate predicted from previously derived effective viscosities. 


Variational Principle for Nonlinear Spiral Density Waves 

Shu, Yuan, and Lissauer (1985a) and Borderies, Goldreich, and Tremaine (1985) 
derived a dispersion relation and a wave amplitude relation for nonlinear spiral density 
waves. We have found a Lagrangian functional that contains the same information as 
their averaged equations of motion (Stewart 1986). Both the nonlinear dispersion rela- 
tion and the wave action transport equation can be easily obtained from the Lagrangian 
by the variational principles described in Whitham’s (1974) book. Our formulation 
greatly simplifies the derivation of these relations and goes a long way toward clarifying 
the relationship of nonlinear spiral density waves to other nonlinear wave phenomena 
described in the general literature. We have added an isothermal pressure term to the 
Lagrangian to facilitate comparisons with earlier theories of linear spiral density waves 
in spiral galaxies. 

In addition to the pedagogical value of this work, our Lagrangian formulation pro- 
vides the means of systematic improvements of the published theory. An equivalent 
Hamiltonian field theory can be derived from the Lagrangian. In this formulation it 
becomes clear that Shu et al.’s (1985a) averaging procedure is equivalent to the standard 
mean field approximation that is widely used in solid state physics problems. Systematic 
improvements to the mean field theory that are commonly used in solid state physics may 
be adapted to the spiral density wave problem using our formalism. Since waves in 
planetary rings are dissipated by inelastic particle collisions, it is desirable to add dissipa- 
tion to the Lagrangian formulation. We are currently investigating the possibility of 
adding dissipation by combining the Lagrangian formalism with a Fokker-Planck model 
for diffusion of the radial displacements and velocities about the mean motion in the 
wave. A Fokker-Planck approximation is feasible for this problem because the relative 
velocity between colliding particles is much less than the systematic velocity in nonlinear 
density waves. Since Fokker-Planck models are commonly employed to model wave 
dissipation in condensed matter physics, we are attempting to borrow the refined tech- 
niques that have been developed in that field. Our goal is to formulate a useful 


alternative to the techniques used by Shu et al. (1985b) and by Borderies et al. (1986) to 
model collisional damping of nonlinear density waves in planetary rings. 


IB. ANTICIPATED CAPABILITIES AND GOALS FOR 7/1/87 - 12/31/87 


Due to the unique capabilities of our proposed numerical simulation procedure 
currently under development, we expect to be able to perform meaningful tests of recent 
analytical formulations of particle disk dynamics. Preliminary results for an unperturbed 
disk have already been obtained: Computed values for the pressure tensor components 
are in good agreement with analytical formulae derived in Shu and Stewart (1985). 

By simulating a wide range of particle number densities and perturber mass ratios, 
we hope to determine general criteria for moonlets to clear gaps or to produce observable 
disturbances in planetary rings. Previously derived criteria were based on very primative 
models of the collisional damping process and limited by the assumption of complete 
damping of forced eccentricities between successive encounters with the perturber (eg. 
Greenberg 1983). 

Given the detail of the observations of the wavy edges in Encke’s division of 
Saturn’s rings (Cuzzi and Scargle 1985), we hope to make special efforts to produce 
detailed models of moonlet-perturbed edges of planetary rings. The most recent model- 
ing effort (Showalter et al. 1986) completely neglects the effects of particle-particle colli- 
sions that we are attempting to simulate. 

An improved theory for radial mass transport in planetary rings is to be derived 
which avoids the assumptions of previous derivations. In particular, there is no necessity 
to identify the off-diagonal component of the velocity dispersion tensor with the effective 
viscosity. 

The propagation of nonlinear spiral density waves in rings is to be described by a 
Lagrangian field theory. Higher order corrections to the approximate theory of Shu et al 
(1985a, b) and Borderies et al (1985, 1986) can be derived with this formulation. An 
alternative method of modeling collisional damping is to be formulated which is better 
able to handle large azimuthal wave numbers characteristic of moonlet perturbed sys- 
tems. 

We hope to place strong constraints on the degree of runaway growth during plane- 
tary accretion. This work may lead to the first simulation of the earliest stages of plane- 
tary accretion which avoids the isotropic "particle-in-a-box" approximation used by 
Greenberg et al. (1978) and others that have been criticized by Wetherill and Cox (1984, 
1985). 
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Abstract 

A self-consistent set of equations for the velocity evolution of a general plan- 
etesimal population is presented. The equations are given in a form convenient for 
calculations of the early stages of planetary accumulation, when it is necessary to 
model the planetesimal swarm by the methods of gas dynamics, rather than follow 
the orbital evolution of individual bodies. To illustrate the relative importance of 
the various terms of these equations, steady state velocities of a simple planetesimal 
population, consisting of two different sizes of bodies, are calculated. Dynamical 
friction is found to be an important mechanism for transferring kinetic energy from 
the larger planetesimals to the smaller ones, providing an energy source for the small 
planetesimals that is comparable to that provided by the viscous stirring process. 
When small planetesimals are relatively abundant, gas drag and inelastic collisions 
among the smaller bodies are of comparable importance for dissipating energy from 
the population. 
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